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Abstract. The body and spatial representations of rigid body motion corre- 
spond, respectively, to the convective and spatial representations of continuum 
dynamics. With a view to developing a unified computational approach for 
both types of problems, the discrete Clebsch approach of Cotter and Holm J^) 
for continuum mechanics is applied to derive (i) body and spatial representa- 
tions of discrete time models of various rigid body motions and (ii) the discrete 
momentum maps associated with symmetry reduction for these motions. For 
these problems, this paper shows that the discrete Clebsch approach yields 
a known class of explicit variational integrators, called discrete Moser-Veselov 
(DMV) integrators. The spatial representation of DMV integrators arc Poisson 
with respect to a Lie-Poisson bracket for the semi-direct product Lie algebra. 
Numerical results are presented which confirm the conservative properties and 
accuracy of the numerical solutions. 
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1. Introduction 

The Hamiltonian structure of continuum mechanics in the material, inverse ma- 
terial, spatial and convective representations was introduced in Holm, Marsden and 
Ratiu 1|12) . This work identified the body and spatial representations of rigid body 
motions as prototypes for the respective convective and spatial representations of 
continuum dynamics. Its comparison of the spatial and convective representations 
also put the Hamiltonian treatments of elasticity by Holm and Kupershmidt l|l0l) 
and by Marsden, Ratiu and Weinstein l|l9! ) into a unified framework. 



The convective representation. The convective and also the inverse material (aug- 
mented Eulerian) representations offer alternative descriptions of continuum mod- 
els. The motivation for the convective representation of the continuum comes arose 
from a number of sources in the 1980s, including the study of relativistic adiabatic 
fluids by Holm ifllh . stability analysis of the coupled rigid body-beam and plate 
models of Krishnaprasad and Marsden iflih and the geometrically exact rod and 
plate models of Krishnaprasad, Marsden and Simo 1)161) . 

Semi-direct products. Holm, Marsden and Ratiu l)l3|) derived the Euler-Poincare 
(EP) formulation of the Eulerian fluid equations for an ideal fluid by applying sym- 
plectic reduction to Hamilton's principle for fluids. Legendre-transforming the EP 
theory recovered the semidirect-product Lie-Poisson Hamiltonian theory that had 
been discovered and applied earlier for nonlinear stability analysis in Holm, Mars- 
den, Ratiu and Weinstein l[l4j) . A key step in the analysis of nonlinear stability 
of fluid equilibria relics on the existence of "Casimirs" - quantities whose Lie- 
Poisson bracket vanishes with all Eulerian (spatial) fluid variables because of right - 
invariancc of the Eulerian variables under rcparametcrisation of the Lagrangian 
labels. Because their Poisson brackets with the Hamiltonian vanish, the Casimirs 
are conserved quantities. We shall use this result to verify that our numerical ex- 
periments preserve the Lie-Poisson structure for the problems we consider below 
by explicitly showing that the values of the corresponding Casimirs are preserved. 



Circulation theorems. Fluid mechanics literature widely refers to the rcparametcri- 
sation of labels as fluid parcel relabelling and attributes the existence of the Kelvin 
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circulation theorem for ideal flow to the application of Noether's theorem for the 
particle relabelling symmetry group. Holm, Marsden and Ratiu l)l3|) showed that 
when advected quantities are present, a corollary of the EP framework is a geo- 
metric form of the Kelvin circulation theorem referred to as the Kelvin Nocther 
theorem. In this framework, Holm, Marsden and Ratiu l(l^)l(l3j) further revealed 
the utility of simple finite dimensional examples, such as the heavy top, by demon- 
strating that they also exhibit a Kelvin Noether theorem. This theorem together 
with the EP equations form an essential ingredient in geometric models of idealised 
continua. 

Variational integrators. Geometric numerical methods seek to transfer these pow- 
erful concepts in geometric mechanics to computational models. In 1991, the pi- 
oneering work of Moser and Veselov |2fih revealed integrable classical mechanical 
systems which have integrable discrete time counterparts. They considered the 
free rigid body as one example and derived a discrete analogue to the Euler- Arnold 
equations for rigid body motion in the body description. These integrators, referred 
to as discrete Moser- Veselov (DMV) integrators, conserve the rigid body energy to 
an arbitrary order of the time step size and angular momentum to numerical round 
off. 

Moscr's and Veselov's key step was to form a discrete Hamilton's action principle 
and then derive "variational integrators" which preserve the variational structure. 
Although the number of contributions to this approach is too extensive to list here, 
the reader may follow some important aspects of its development in Bobenko and 
Suris 10), Marsden, Pekarsky, and Shkoller . Marsden and West (H2), McLachlan 
and Scovel |24h, W endlandt and Marsden l(3fl ) and, more recently, Leok, Marsden 
and Weinstein l)l7|) who provide a differential geometric foundation for variational 
integrators applied to mechanical systems. The number of numerical studies sup- 
porting the theory appears less extensive, however, largely due to the absence of a 
unified computational framework for deriving practical variational integrators. 

Practical integrators. DMV integrators are explicit. Cardoso and Leite (0) cast 
the expression for the discrete angular momentum of Moser's and Veselov's rigid 
body into an algebraic Ricatti equation and solved it by Schur decomposing the 
Hamiltonian matrix. This gives a nearly explicit DMV algorithm, except for the 
iterative Schur decomposition. McLachlan and Zanna (|25h recently demonstrated 
how to avoid the costly computation of the Schur form by using an explicit spectral 
decomposition of the Hamiltonian instead. The Hamiltonian can be decomposed 
in this way whenever its characteristic polynomial can be solved analytically. The 
simple models considered here did not require this optimisation step. 

Aims. In this paper, we seek to develop a unified practical computational frame- 
work for continuum dynamics by deriving a class of explicit variational integrators 
for various rigid body motions in the body and spatial representations using the 
same underlying approach, referred to as the discrete Clebsch approach. The con- 
tinuous time Clebsch approach provides a systematic means of deriving the EP 
equations from a symmetry reduced Hamilton's action principle and the "momen- 
tum maps" (see (|2 lr i for an explanation of this term) associated with this symmetry 
reduction. For the case of rigid body motions, we aim to show that the discrete 



3 



Clebsch approach provides discrete analogues, the discrete EP equations and the 
corresponding momentum maps. We also aim to assess the conservative properties 
and accuracy of the DMV integrators derived in this framework through numerical 
experiments presented herein. 

1.1. Approach. We begin in section El by returning to arguably the most cele- 
brated model taught in mechanics, the classical free rigid body in continuous time 
whose description, in the body frame, can be found in an exceptionally lucid and 
concise form in Marsdcn and Ratiu's introductory text book on geometric mechanics 
|2lL We slightly modify the notation for the purposes of describing the kinematics, 
symmetries and equations of motion of the free rigid body in discrete time. We 
also include the spatial description of the rigid body. In each case, we apply the 
Clebsch approach of Cotter and Holm (0) , a discrete extension of the earlier work 
of Cendra and Marsden 10) and Holm and Kupershmidt iflpfo . to both derive the 
basic discrete EP equations in body variables and the discrete EP equations with 
an advected parameter in spatial variables. 

We also derive the right and left infinitesimal equivariant momentum maps as- 
sociated with cotangent lifted left and right actions of SO(3) on the canonical and 
augmented cotangent bundle respectively. We finally demonstrate the conservative 
properties and accuracy of the discrete EP equations through numerical experi- 
ments of free rigid bodies, heavy tops and coupled rigid bodies in section [5] 

Appendix A provides a summary of the main features of the spatial and body 
representations of each of the models that we consider, together with a comparison 
of their continuous time counterparts. The geometric form of the discrete and 
continuous time models, derived from a Clebsch approach, are remarkably similar. 
For completeness, appendix B provides the spatial version of McLachlan's and 
Zanna's (unoptimised) DMV algorithm l|25h . We now state the main results of this 
paper. 
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1.2. Summary of main results. 

Discrete Clebsch approach. Wc show how this discrete Clebsch approach not only 
recovers the body representation of the basic discrete EP equations of l|26h but also 
yields a spatial representation of the EP equations with an advectcd parameter 
which correspond to the Lie-Poisson equations on the dual space of a semidircct 
product lie algebra, discovered by Bobenko and Suris J3) for a more general class of 
problems when the phase space is G x G. Bobenko and Suris also showed that the 
Lax representation in the spatial frame is gauge invariant to its representation in the 
body frame. Bobenko and Suris's discovery of the semi-direct product in discrete 
time enables the transfer of the significant work of Holm, Marsden and Ratiu on 
semi-direct products in continua to the discrete time case. Prior to their paper, 
the relation between EP equations with an advectcd parameter and Lie-Poisson 
Hamiltonian systems required invertibility of the partial Lcgcndre transform. This 
work overcame this restriction by developing the EP theory entirely within the 
Lagrangian framework. 

Momentum maps. We derive right and left momentum maps for the respective 
left and right reductions of the Lagrangian by SO (3). The right momentum map 
is associated with the cotangent lifted left action of SO (3) on the approximated 
cotangent bundle and its image is the spatial angular momentum, which is an 
element of so(3)*. Left invariance of the action principle implies conservation of 
spatial angular momentum m. The image of the left momentum map is the body 
angular momentum M which is only conserved when the action principle is right 
invariant. This is the limiting case when one or more of the moments of inertia are 
equal and is not considered in this paper. 

Conserved quantities. We prove that both representations of the DMV integrator 
conserve the spatial angular momentum and demonstrate through numerical exper- 
iment that the DMV algorithms conserve spatial angular momentum to numerical 
round off. We also demonstrate the correct computation of the Casimirs 1 1 AiT 1 1 § 
and {||/||2, det(I)} of the Lie-Poisson bracket for the respective body and spatial 
representations and that the numerical solution conserves energy to an order of the 
time step. 

The heavy top. We apply the same approach to the classical model of the heavy top. 
For a detailed treatment of the motion of the heavy top we refer the reader to the 
work of Lewis, Simo and Marsden l(l8t ). They consider three types of heavy tops, 
the asymmetric, the tilted and sleeping Lagrangian tops and extend the classical 
studies of the heavy top by using the reduced energy momentum to derive relative 
equilibria. 

A result which seeks numerical validation is the observation that stable branches 
of steadily precessing Lagrangian tops bifurcate from the branch of sleeping La- 
grange lops throughout the range of angular velocities for which the sleeping top is 
stable. Wc do not address the extension of these results to the discrete Lagrangian 
framework but instead present some geometric properties in the discrete framework 
(in both representations) together with numerical results of the motion of the heavy 
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top in the body representation only. Appendix A provides a description of the geo- 
metric properties of the body and spatial representations of the heavy top motion 
in discrete and continuous time. 

The coupled rigid body. We apply the discrete Clebsch approach to the classical 
coupled rigid body. We begin by presenting a slightly revised description (in the 
frame of body 1) of the M 3 reduced coupled rigid body which is clearly described in 
the thesis of Patrick (j^) and the concise paper by Grossman, Krishnaprasad and 
Marsden (0). The last table of appendix A provides a spatial description of the 
continuous time problem together with a corresponding description in discrete time. 
We also hope that this description of the discrete coupled rigid body will be ex- 
tended to form a discrete analogue of the engineering related work on coupled rigid 
bodies and geometrically exact rods by Simo, Posbergh and Marsden Ij28l) and the 
study of the dynamics and stability of the coupled rigid bodies by Sreenath, Krish- 
naprasad and Marsden l)29{) . We provide some results from numerical experiments 
of the coupled rigid body as viewed in the frame of body 1. 

1.3. Important related works. Similar approaches derive discrete equations of 
rigid body motion for optimal control problems. Bloch, Crouch, Marsden and 
Ratiu (|l|), for example, derive the symmetrised rigid body equations by introducing 
optimality constraints in the action principle. We distinguish our approach from 
theirs in two ways. Firstly, although they consider the rigid body motion as an 
optimal control problem with an associated constrained action principle, they do 
not identify the constraints as Clebsch variables and derive the momentum maps. 
Secondly, they present left and right trivialisations of T*SO(3) where as we present 
body and spatial representations of a left SO (3) action invariant Lagrangian only. 
The authors make this point when distinguishing their approach from that of Holm 
and Kupershmidt l)lO|) . We use the expression for the (left) momentum map to 
prove that the flow on the cotangent bundle preserves spatial angular momentum 
and derive the equations of motion. 

The subsequent presentation of the discrete equations of motion and momentum 
maps for the heavy top appears unique in so far as the only citation found on 
geometric integrators for the heavy top, by Celledoni and Safstrom (Q), does not 
consider DMV integrators nor the spatial representation. Work on practical aspects 
of geometric integrators for the discrete coupled rigid body problem is not cited in 
the literature although Patrick references several of his own Maple scripts for 
computing various properties of the coupled rigid body. 

2. The Free Rigid Body 

In this section, we slightly modify the description of the free rigid body given 
in chapter 15 of (|2lT ) as a discrete time problem in both the body and the spatial 
representations. Consider a free rigid body as a solid body occupying a reference 
configuration B S R which is free to move in the container C = R 3 by rotations 
about a fixed point. Material points t S B are position vectors whose components, 
relative to a fixed orthonormal basis (Ei, E2, E3) in B, are the material coordinates. 
A configuration € of B is a continuous, invertible and orientation preserving map 
4> : B — > C from material points to spatial points in the container. The spatial 
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points are position vectors whose components, relative to (ei, 62,63), the right- 
handed orthonormal basis of C, are spatial coordinates. 

Discrete time dependent families of configurations of B are referred to as a dis- 
crete time motion of B and are written in terms of the forward map x = "ip{£, tk) = 
ipt k (£), k G Z + , which in addition to having the properties of the configuration also 
satisfies ijj(£,t ) = ipt (fy = 

This last property together with rigidity of the body and continuity of the motion 
implies that the configuration of B may be identified with SO (3) and the forward 
map is written as 

x = ipk{£) = A^, A fc G 5*0(3), (1) 

where, for notational convenience, tpk '■= ipt k and := A(t/ £ ), represents the 
attitude of the body. Put simply, the forward map is the position of a label £ in 
the container. The value of the label is its position in the container at time to. 

The body coordinates of a material position vector are its components relative to 
a time- dependent basis (£1, £2, £3 )(£*:) which is defined by £i(tfc) = AkE^, i := 1 — * 3 
and hence is attached to the rigid body that rotates about the origin of C. 

2.1. Discrete velocities. The image of the backward or inverse map I = ip^ (pc) 
is the label whose position at time to = is x. With the configuration identified as 
50 (3), this is given by 

£ = ^7 1 (x)=A- 1 x. (2) 

It follows that label positions at consecutive times tk-i and tk are related by the 
composite of forward and backward maps which are group products 

M*) = fa V'fe-i 1>h-i(t) = AfcA^Afc-iW- (3) 

The product of a matrix and an inverse matrix at consecutive times is referred 
to as a discrete velocity. Just as the (continuous time) spatial and body velocities 
are right and left invariant respectively, so too are their discrete counterparts. 

The spatial discrete angular velocity uJk+i and body discrete angular velocity 
fife+i respectively satisfy the reconstruction formulae 

Afe+i = ujk+iAk, 
A/c+i = Afcfifc+i. 

The two discrete velocities at time tk+i are related by 



tt k+ i = A^Wfe+iAfc = Ad A Tw k +i. (5) 

Remark 2.1.1. The reconstruction formulae conserve labels. This statement is 
confirmed by considering the spatial reconstruction formula 



Al +1 



Xk+l ~ LUk+lXk 

ATfc+i = h^Xk 

£k+l = t-h- 



(G) 



With the spatial and body discrete velocities defined, we now consider the geo- 
metric mechanics of rigid body motion in discrete time. 
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3. Discrete Constrained Variational Principle 

Moser and Veselov |2^l consider a Lagrangian L:GxG^K which is a smooth 
map defined as 

L(A fc , A fc+1 ) = Tr(A k I Al +1 ) - Tr (9 fe+1 (A fc+1 A^ +1 - I d )) . (7) 

The constant inertia matrix Iq G V* — 52 (M 3 ) is a matrix of the positive symmet- 
ric covariant two-tensors on R 3 . We identify V with V* using the metric Tr(A T B) 
so that V = 5 2 (M 3 ) is the 3x3 matrix of symmetric covariant two-tensors on R 3 
dual to 52 (R 3 ). The pairings in the above Lagrangian are therefore between V and 
V*. 

We denote the projection of elements g G G onto 5 2 (R 3 ) as sym(g), where symQ 
is a projection operator defined as sym(g) := |(<? + g T ). The symmetric Lagrange 
matrix multiplier Ofe+i G V* enforces orthogonality of Afc 1 . 

The problem is to find the sequence {A;., Afc + i} which satisfies the discrete sta- 
tionary action principle 2 

= 5S d = J2 L (^k,Ak+i)- (8) 

k 

A necessary condition for extremising this functional is that {A^, Afc + i} satisfy the 
discrete Euler-Lagrangc equations. We shall now derive the solution in the spatial 
and body representations. 

4. Symmetry Reduction 

Body representation. The discrete differential geometric formulation of the DMV 
system is given in lfl7L The principal G-Bundle (G x G, G, n) together with the nat- 
ural projection n: GxG^GxG/G furnish the description of discrete symmetry 
reduction to the body representation. 

Definition 4.0.2. The (left) diagonal action of G on G x G is defined as ^ : 
Gx(GxG)^GxG\ *(/, ( 5 , h)) = f ■ (g, h) = (fg, fh). 

The discrete Lagrangian in equation is invariant under the (left) action of \&. 
As depicted by the two curved arrows in figure 1, one may reduce the Lagrangian 
on the G-bundle by this action to obtain the reduced Lagrangian L : G — > R given 
in body variables by 

L(fijfc+i) = Tr(sym(n k+1 )I ) - Tr (6 fc+ i(fi fe+1 ^ +1 - I d )) . (9) 



This constraint is required to derive the basic discrete EP equations on SO (3) and not on 
S 2 (R 3 ). 

2 Holm, Marsden and Ratiu il3f) point out that this is not strictly a variational principle 
(because the variations of the dynamical variables are constrained) but a Lagrange D'Alembert 
principle. 
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ff-^Afc) 




Figure 1. The principal G-Bundle upon which the discrete Lagrangian is de- 
fined. The two curved arrows represent the diagonal action of A k on (A^ , Afc+i). 
The vertical arrow represents the fiber over . 



Spatial representation. To reduce to spatial variables, we must firstly modify the 
definition of the Lagrangian for reasons which are give in section 6. We consider 
the Lagrangian defined on the augmented space, 

L : SO(3) x 50(3) x V* -> M, (10) 

given by 

L(A k ,A k+1 ,I ) = Tr (AfcloAjki) 1 - Tr (e fe+1 (A fe+1 A^ +1 - I d )) . (11) 

We must define how SO(3) acts on V*. Let $*(g)-I, I € V*, geG = 50(3) be 
the left group anti-representation on V* . This action defines an orbit of the inertia 
matrix in V* , which is a dynamical variable in the spatial frame. 

The left action <!>*(#) on V* is given by 

:G xV* -^V* \ <f>*(g) ■ I = gig- 1 . (12) 
The corresponding right action is given by 

: V* x G V* \ I ■ $*(g) = g~ x Ig, (13) 
which is equal to the left anti-representation $*(g _1 ) • /. 

Definition 4.0.3. The (right) augmented diagonal action of G on G x G x V* 

is defined as : (G x G xF) x G -> G x G xV* | V'((g, h, a), /) = (g, h,a)-f = 
(gf,hf,a-Z>*(f)). 
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This discrete Lagrangian, which is now a function of Iq, is invariant under the 
(right) action of g G 50(3) 

L(A k g,A k+1 g,g T I g) = L(A k , A k+1 , I ). (14) 

Reduction of the Lagrangian by the augmented diagonal action gives the reduced 
Lagrangian L : G x V* — > M given in spatial variables as 

L(u>k+i,h) = Tr(sym(uk+i)h) - Tr (6fc + i(wfc+iwf +1 - Id)) ■ (15) 

5. Clebsch Potentials and Momentum Maps 

The symmetry reductions have an associated momentum map which may be 
derived by adopting the discrete Clebsch approach of Cotter and Holm (0). An 
associated or augmented Lagrangian is defined by adding a Clebsch potential to 
enforce the reconstruction formula. We will consider the approach in body and 
spatial representations. 

5.1. The body representation. In body variables the augmented Lagrangian is 
given by 

Tr 

L' := Tr (I oSy m(n k+1 ))- — (P fc T +1 (A fe+1 - A k n k+1 ))-Tr (e fe+1 (O fe+1 fi£ +1 - I d )) . 

(16) 

The stationary constrained action principle gives Clebsch relations which include 
the discrete Euler-Lagrange equation 

V Afc £'(A fc _i, A fe ) + V Ak L'(A k , Afc+i) = 0. (17) 

Evaluation of the derivatives of L' and subsequent rearrangement gives an evolution 
equation for the Lagrange multiplier P k 

Pk+i = Wfc+i-fV (18) 
This equation together with the discrete evolution equation for A^ describe the 
discrete flow on the cotangent bundle and are referred to as the discrete symmetrised 
equations and are analogous to the continuous symmetrised equations. 

The Clebsch relation paired with the variation SQ k +i in the discrete action prin- 
ciple gives the expression 

J* +1 = skcw(A^P,), (19) 

which satisfies the definition of the momentum map for cotangent lifted left actions 
of SO(3) 



= (P fe oA fc ,C), 

where the bilinear operation o : V x V* — > g* is defined in (H3) and (A, B) = 
-\Tr{A T B) for any matrices A,B eV. 

It is well known that momentum maps for cotangent lifted actions are equivariant 
and consequently Poisson. The image of this map is the body angular momentum 
Affc+i := ^o^fc + i — ^fc+i^o- 
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Substitution of the symmetrised equations for the discrete flow on the cotangent 
bundle 



P k+1 = P k Q k+1 , 

A -AO ( 21 ) 

into the right momentum map gives the discrete basic EP equation defined on 
50(3) for rigid body motion in the body frame 

M k+1 = Ad* n -rM k . (22) 

5.2. The spatial representation. We now derive the equations of motion in the 
spatial representation by applying the discrete Clebsch approach. In spatial vari- 
ables, the Clebsch modified discrete Lagrangian is given by 

Tr 

V := Tr (I k sym(u k+ i)) - — (P fc T +1 (A fc+ i - uj k+ iA k ) 

Tr (23) 
- —{Jk+i{h+i - w fc+ i4Wfc +1 )) - Tr (Q k+1 (uj k+1 ujl +1 - I d )) . 

The stationary constrained action principle gives Clebsch relations which are the 
discrete Euler-Lagrange equations 

VA k i , (Afc_i,A fc ) + VA k i , (Afc,A fc+1 ) = 0, / 

(24) 

Vi k L'(I k _ 1 ,I k ) + V Ik L'(I k ,I k+1 ) = 0. 

Evaluation of the derivatives of L 1 and subsequent rearrangement gives, respec- 
tively, 

Pfc+1 = UJ k+1 P k , 

Jfc+i = uj k+ i (-sym(u k ) + J k )uj k+1 . (25) 

G k 

Additionally, the Clebsch relation paired with 5oj k +\ gives 

I k+1 uj k+1 + A k P k+1 LL! k+1 + I k uj k+1 J k+ iUJ k+ i = 9fe+i. (26) 

Using the symmetry property of <d k , the equations P? = P£ +1 uj k+1 and G k = 
u> k+1 J k +iLi! k +i gives an expression for J k+1 

Jt +1 =skcw(P k Al) + [G k7 I k }, (27) 

which satisfies the definition of the momentum map for cotangent lifted left actions 
of S0(3) 

(Jt+iX) = <ft,£ C A fc ) + (G k ,2 c I k ) 
= (P k oA k + G k oI k ,Q. 

The image of J^ +1 is the spatial angular momentum m k+ i := I k+ iuj k+ i — 
The spatial representation of the discrete EP equations with an ad- 
vected parameter are 

li 



TOfc+i = Ad* uT m k -Vi k LoI k , (29) 

4+i = w k+1 I k ul +1 . (30) 

Lemma 5.2.1. The spatial angular momentum is a conserved quantity. 

Proof. Substituting the symmetrised equations for P k and into the expression 
for the left momentum map gives 

?7i fc+ i = WkiPk-iAl^ - A k -iP k _ x )oj k - [h, G k ] 

= w k m k u k + uJ k [I k -i,Gk-i]u k - [Ik,G k ] 

Using G k = LOkGk-iUjl - sym(uj k ) and I k = Lj k I k -i^ k gives 



(31) 



m k+1 = io k m k tJ{ + [h,Uk + 

= U} k m k ul + I k u k - uj k lk - Lo k h + I k 0J k 
= io k m k ujJ. +m k - LU k (I k uJk - LJ k I k )u k 
= m k . 



(32) 



□ 



So for conservation of spatial angular momentum, the coAdjoint orbits of the 
action of SO(3) on so(3)* must take the form 

Ad^Tino = m + V/ & £ o I k . (33) 

6. Poisson Brackets and Semidirect Products 

Bobenko and Suris (0) show that the right reduced discrete EP equations are 
Lie-Poisson w.r.t. a semidirect product lie algebra. We show how the spatial 
representation of the discrete EP equations for the rigid body, given above, are 
related to Bobenko's and Suris's result, which is generalised to a class of systems in 
which the heavy top is one example. We will firstly present the Lie-Poisson bracket 
for the rotating rigid body in the spatial frame by excluding the centre of mass 
vector x from the study of the heavy top in the spatial frame performed by Holm, 
Marsden and Ratiu Ijl2h . 

Unless otherwise stated, we omit the time index subscripts for ease of notation. 
We point out a few minor differences in our notation with Bobenko's and Suris's. 
The first point is that Bobenko and Suris consider a dynamical variable p £ V, in 
the linear space in which the group is represented, whereas we consider a dynamical 
variable I G V* , m the dual of this linear space. Our notation is consistent with the 
work by Holm, Marsden and Ratiu on semidirect products in the EP description of 
the continuum l|13|) . 

We extend the definitions given in section 4. The corresponding left anti- 
representation of the lie algebra on V* is 

<T : 3 xV*^V* | $*(C) -I=[C,I], (e 9 = so(3), (34) 

with the right anti-representation following from the definition of the right group 
action on V* 



: V* x g V* |/.**(C)J=[/,C], Ceg, 
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(35) 



where it follows that 

$*(CW = -**(C T W = -^-$*(0- (36) 

Unlike the material representation of the general class of Lagrangians on G x G 
that Bobenko and Suris consider, for any group G, the material representation of 
the Lagrangian for the rigid body on 50(3) x 50(3) is not invariant under the 
diagonal right action of G = 5*0(3). The only exception to this occurs when one 
or more of the principal moments of inertia are equal to each other for which the 
isotropy group G/ for I is the subgroup SO(2, 1) or 50(3) respectively. We do 
not consider this special case and instead consider the Lagrangian defined on the 
augmented space, given in equation 1111 which is invariant under the augmented 
diagonal right action defined in section 4. 

The discrete EP equation with an advectcd parameter 1^ given in cquation l29l is 
defined on G x V* . Under the assumption of invertibility of the partial "Legcndrc 
transformation" 

(wk+i.ifc) G G x V* -> (mfc+i, J fc+ i) e g* x V\ (37) 

given by Bobenko and Suris (equation 4.19, 10)), equation defines the smooth 
map 

(m k ,I k ) ^ (m fc+ i,4+i) EQ*xV*. (38) 
This map is Poisson w.r.t. to the (±) Lie-Poisson brackets on the dual of the 
semi-direct product algebra [s* = 0*(S)U*]± which have the form 

{F 1 ,F 2 } ± (m,I) = -^Tr{±m T [\7 m F u \7 m F 2 ] 

T I (V/F 2 ■ $(V mJ F\) - y I F 1 ■ $(V m F 2 ))}, 

where $ is the right representation of q in V which is related to the anti-representation 
by the expression given in (equation 3.16, (0)). These are the (±) Lie-Poisson brack- 
ets on s* for the right representation $ of g on V given by Holm, Marsden and Ratiu 
(section 5.5, lfl3) ) for the spatial representation of the heavy top (in the absence of 
a gravitational potential). The same authors also give these brackets for the more 
general case when G is any group which acts on V* from the right in (equation 
2.14, ©). 

The (-) Lie-Poisson bracket in equation[2Hlis the same bracket defined by Bobenko 
and Suris in (equation 4.20, (Q)) and so we conjecture that the proof (page 12, (Q)) 
of the Poisson property of the map given in equation 1381 applies here. Note that 
they state the (-) Lie-Poisson bracket of the semi-direct product lie algebra cor- 
responding to minus the left anti-representation of G on V* , whereas we state 
the form of the (±) brackets corresponding to the right representation of 5*0(3) 
on V. The relationship between the representations given in equation 1361 permits 
the simple interchange between the bracket corresponding to the (-) left and right 
representations. The right representation is a prototype for idealised fluids. 

We close this section by commenting on the Casimirs of this Lie-Poisson bracket 
which are less well-known than for the standard rigid body bracket. We firstly 
express the Lie-Poisson bracket in a more concise form 

{F 1 ,F 2 } ± {m,I) = ±(m, [V m Fi, V m F 2 ]> ± (I o V,F 2 , V m Fi> =F (Jo Vj.Fi, V m F 2 ). 

(40) 

Holm, Marsden and Ratiu l|12l ) show that the Casimir functions on the Poisson 
manifold [s*]± are the set of functions invariant under the coAdjoint action of the 
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Lie group S = SO(3)(§)V* . For the right representation of A on V* , the definition 
of this action, given in a general form in (equation 2.09, is 



(A, J)(C, I) = (Ad* AT ( + J ■ $(A T ) o / • $(A T ), I ■ $(A T )}. 



(41) 



It follows that invariants of / under conjugation by A, such as det(I) and ||/||2 
are invariant under this coAdjoint action and are subsequently Casimirs of the 
Lie-Poisson bracket. 



In this section, we extend the main results of the spatial and body representations 
of the rigid body, presented thus far, to the heavy top. We consider the kinematics 
and symmetries of the tilted Lagrange top, described in the body frame by the 
orientation T k of the vertical axis z and the body angular velocity (See figure (2)) 
and in the spatial frame by the position of center of mass \k relative to the support, 
the inertia tensor and the spatial angular velocity (see figure (3)). 

The Lagrangian top is a special case of the heavy top which has a symmetric 
inertia matrix and its centre of mass lies on its axis of symmetry. The motion of the 
tilted Lagrangian top is distinct from the sleeping top in that it both spins about 
its axis of symmetry and precesses about each of the spatial axes, maintaining a 
positive vertical coordinate. The axis of symmetry of the sleeping top remains 
parallel to the gravitational vector — gz, however. 

7.1. The body representation. 

It is well known that the motion of the heavy top breaks the symmetry of the La- 
grangian when precession modifies the gravitational potential of the top or its spin 
is not purely about its centreline. The symmetry group for the heavy (Lagrangian) 
top is thus S 1 x S 1 . From Noether's theorem we may deduce that spatial angular 
momentum will only be conserved if the group action for the motion is that of this 
symmetry group. 

The Clcbsch variable constrained Lagrangian for the discrete time heavy top 
motion in the body representation is 



- (J fc+ i,r fe+1 - n k+1 r k ) + Tr (e k+1 (n k+1 nl +1 - i d )) , 

where h is the fixed time interval. The (right) momentum map corresponding to 
the left augmented diagonal action of 50(3) on 50(3) x SO(3) x R 3 is 



7. The Heavy Top 




(42) 



fc+i — 



PfeoAfc + r fc oj fe , 



(43) 



with 3 k = J fc + mgxo- 

The Clebsch relations include the discrete Euler-Lagrange equations 



V Afc £'(A fe _i, A fe ) + \7 Ak L'(A k , A fe+ i) = 0, 
Vr fc i'(r fc _i,r fc ) + Vr k L'{r k ,T k+1 ) = 0. 



(44) 
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Figure 2. The heavy top as viewed in the body frame. The heavy body is 
attached to the spatial frame at an arbitrary point (in this diagram this point is 
the origin). The motion is composed of two components, precession and spinning. 
The unit vector z , representing the direction of gravity, rotates about each axis 
of the heavy top with body angular velocity O. Spatial angular momentum is only 
conserved for motions purely about E3, however. The body frame also spins about 
its centreline axis but this is only observable in the spatial frame. The vector \o 
from the point of support to the centre of mass of the top (c.o.m.), which lies on 
the centreline of top, remains fixed in the body frame. 



Evaluation of the derivatives of L' and subsequent rearrangement gives, respec- 
tively, 

1 

P fc+ i =Pfcft fc+ i, 

(45) 

Substitution of these two equations into the equation for the (right) momentum 
map gives the equations of motion in the body variables 

M k+1 = Ad* T M k + rngT k o X 0, 

(46) 

Tfc+i = f2 fc+1 rfc, 

with AI k := j^skew ((V f2(i+1 L) T f] fe+1 ). 
7.2. The spatial representation. 
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Figure 3. The heavy top as viewed in the spatial frame. The heavy body is 
attached to the spatial frame at an arbitrary point (in this diagram this point is the 
origin). The motion is composed of two components, precession and spinning. The 
vector x from the point of support to the centre of mass of the top (c.o.m.), which 
lies on the centreline of the top, rotates about each of the spatial axes with spatial 
angular velocity oj, maintaining a positive vertical component. Spatial angular 
momentum is only conserved for motions purely about 63, however. The body 
frame also spins about its own centreline axis and hence invariance of the Inertia 
matrix under the action of S 1 is also required for spatial angular momentum con- 
servation. 



The Clebsch variable constrained Lagrangian for the discrete heavy top motion 
in the spatial representation is 



Tr 

- -y (Jfc+i(Jfe+i - c*>fe+iifcWfc+i)) - <J£ +1) Xfc+i -LJk+iXk) 
+ Tr (9 fc+ i(c^ + i^+i ~ J d)) ■ 



(47) 



The (left) momentum map corresponding to the augmented right diagonal action 
of 50(3) on SO(3) x SO{3) xV* xM 3 is 



= A k o P k + G k oI k +3*0 Xk, 



(48) 



with J£ := J£ + mgz and G k := u k+1 J k+ iuj k+ i. 
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The image of J^ +1 is the spatial angular momentum which is not conserved 
unless the motion of the top is about the axis parallel to the gravity vector. The 
Clebsch relations give the discrete Euler-Lagrangc equations 

VA fc £'(A fc _i, A fc ) + V Afc £'(A fc! Afc+i) = 0, 

V Ik L'(I k ^,I k ) + \7 Ik L'(I k ,I k+1 ) = 0, (49) 

^Xk L '(Xk-i,Xk) + V Xk L'(xk,Xk+i) = 0. 

Evaluation of the derivatives of L' and subsequent rearrangement gives, respec- 
tively, 

Pk+i = ^k+iPk, 

Jfc+1 = UJ k+1 (-(V Ik L) T + Jfc)^fc+1, (50) 

J fc+1 - ^fc+l J fc- 

Substitution of these three equations into the equation for the (left) momentum 
map gives the equations of motion in the spatial variables 

?7ifc + i = Ad^T-mk - V/ fc L o I k + mgxk o z, 

Ik+l = ^k+lh'-^k+i, (51) 
Xfe+1 = ^k+lXk, 

with m k := -f^skew ((V Wfc+1 L) T u;fc + i). 
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8. The Coupled Rigid Body 




Ball and socket joint 



Figure 4. The coupled rigid body as viewed in the frame of body 1. Each 
body is attached from its centre of mass (c.o.m.) to a ball and socket joint. In 
the IR 3 reduced configuration space, the origin of the spatial frame is the centre of 
mass (C.O.M.) of the coupled rigid body. In the frame of body 1, the motion is 
composed of two components, precession and spinning. The vector Ad^, represent- 
ing the position of the centre of mass of body 2 in the frame of body 1, rotates 
about the origin with relative orientation A — A 1 A2 . (ft ari d V denote the angles 
between the body axes E± and the vertical and 6 denotes the angle between the 
body attachments at the joint. Each body also spins about its axes, but only the 
spin of body 2 is observable. 

We now state the Clebsch variable constrained action principle, derive the mo- 
mentum maps and the equations describing the discrete time free motion of two 
rigid bodies coupled by a ball and socket joint. We choose the origin of the con- 
tainer to be at the position of center of mass of the coupled body (C.O.M.) (as 
shown in figure 4) and let the configuration $ be the attitude of body 1 and 2, 
Ai, A2 <E 50(3) relative to a reference configuration. The basic configuration space, 
under the assumption that the centre of mass of the coupled body is stationary, is 
£ = 50(3) x 50(3). 1 

We denote the total mass as m = mi + m2, the position of the center of mass of 
each body as di and the attitude of body 2 in the frame of body 1, referred to as 
the relative orientation matrix, as A = A^A2. 

It is also useful to define the following terms e = milTt2 , Dij = di® dj , and the 

2 

modified inertia matrix Z; = 7, — ^-Du. Both D12 and Ii are fixed in the body frame 
whereas A is an advected quantity. The discrete Lagrangian on (50(3) x 50(3)) 2 
is 

2 

L = Y,Tr (A,fi 2 Af +lT ) + eTr ((A^ +1 - Aj)L> 12 (A^+ 1 - A*) T ) . (52) 
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8.1. The body representation. The reduced Lagrangian with Clebsch potentials 
is 

2 

L ' = E Tr (liSym(n>Z +1 )) + eTr ((I d - n k+lT )D 12 (I d - Cl k+1 )A k + lT ) 

i=l 

_ Tl (^+i T ( A ^ +1 - A k fl k+1 )) - ^ (j fe+lT (A fc+1 - ^ +lT A fe ^ 2 fc+1 ) A (53) 

-Tr (ef +1 (^ +1 ^+ lT -J d ) 
The Clebsch relations include the discrete Euler-Lagrange equations 

V A? L'(Af- 1 ,A^) + V A? L'(Af,A l ft+1 ) = 0, ie{l,2}, 
V A * L'(A fc_1 , A fe ) + V A * L' (A fc , A fc+1 ) = 0. 
Evaluation of the derivatives of L' and subsequent rearrangement gives, respectively, 



fe+i _ D fc n fc+i 

(55) 

2 

fc+1 



Addition of the Clebsch relations paired with SCI* gives the expression 

2 

J R = J2 skew(Af Pf ) + [J fe , A fe ], (56) 

i=l 

which satisfies the definition of the momentum map for cotangent lifted left actions 
of SO(3) 



1=1 

2 



(57) 



(£2 p i !<>A i + jk oAfe '0- 



i=l 



Lemma 8.1.1. The total spatial angular momentum is conserved by the discrete 
time flow on the augmented bundle. 



Proof. The total spatial angular momentum is 
= m k+1 + m k+1 = ^ AfMf +1 Af 



TO 

i=l 
2 



V 2skew(A k P kT + A k J k Af + k k 2 J kT kf) 

(58) 

2 

^2 S fce W (Af- 1 affif p k ~ lT ) 



m\ + to 2 = TO fc . 



□ 

1!) 



Substituting the expressions in equation [S3] into the right momentum map gives 
the discrete EP equations in body variables 



M k+1 = Ad* nf M k , *G{1,2} (59) 

together with the evolution equation for the relative orientation matrix in the frame 
of body 1 

A fc+1 = Q k+lT A k Q, k+1 . (60) 

8.2. The spatial representation. The Clebsch variable constrained Lagrangian 
for the spatial representation is 



L ' = E Tr (%svm(">t +1 )) + &r [(cj k+1 - I d )D k 2 (u, k+1 I d ) 

i=l 

- Tr (Vf + 1 " (A k+1 - ui k+1 A k )j - Tr ( J k + 1 " (if + 1 - u k+1 I k uj k+ ^ T )) (61) 
k+i T (n k+i _,.,fe+l n fe ,.,k+i T \\ _rr^ fofc+i/ / .,fc+i T ,.,fc+i 



- Tr (J fc+i (Df+ l - u1 +l D1 2 lj« +1 ) 1 - Tr (9f +i K +i - 7 d ) 

where I k = A k IiA kT and D^ 2 = A^D^A^ are now time dependent. 
The Clebsch relations include the discrete Euler-Lagrange equations 

V A *Z'(A*-\ A*) + V A? L'(A*, Af +1 ) = 0, » e {1,2}, 

V^L'iit 1 ,^) + ^ifL'{ll I k+1 ) = 0, i e {1, 2}, (62) 

Vz^L'piT 1 , -°12) + V Df2 L'K 2 , Dfr 1 ) = 0. 

Evaluation of the derivatives of L' and subsequent rearrangement gives, respec- 
tively, 



fc+i _ , ,fc+i pfe 



J* + 1 = u,™ {-sym{u k ) + J k ) w « 

v v ' (63) 



G ■ 



J fc +! = W *+ 1 J fc w£ +1 



Addition of the Clebsch relations paired with 517/ gives the expression 

2 

J L =Y J skew{P k kf) + [G k J k ] + skew([J k ,D k 1 }), (64) 

i=l 

which satisfies the definition of the left momentum map for a cotangent lifted right 
action of SO(3) 



(J L X) = ]T(/f,£ C A.f) + (Gl&c%) + (J k ,^D k 21 ) 

i=l 
2 

= (J2P?otf + G k o% + J k oD k 1 ,0- 
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The total spatial angular momentum m = mi + m-i is the image of J L where 
the spatial angular momentum of each body is 

m\ = 2s/fcew(i 1 X' +1 - e^ +lT (^ +1 - Id)D$), 

pfc+1 

T 1 ( 66 ) 

m\ = 2skew{I^ +1 - eoj% +1 - h)D\ 2 ). 

pfc+1 

1 2 

In summary, the discrete EP equations of motion in the spatial frame are 



i=i 

fk+1 _ , fc+1 fk, ,k+l T 

n k+i _ k+i n k k+i T 



(67) 



We remark that these equations take a very similar form for the n multi-body in 
which the summation is over 1 — » n. Tables 5 and 6 of appendix A summarise the 
main results of this section. We present results from the implementation of the 
body representation of the discrete EP equations with advected parameters in the 
next section. 



9. Numerical Experiments 

This section presents results demonstrating the conservative properties and ac- 
curacy of the rigid body, heavy top and coupled rigid body integrators. The com- 
ponents of the body momentum are compared with the analytic solution for the 
rigid body only, and the Matlab Ode45 integration of the Euler- Arnold ordinary 
differential equations and their variants for the heavy top and coupled rigid body. 
The tolerance of the Ode45 routine is set to le-15. The time step for all numer- 
ical experiments is At = 0.1. Although the figure captions give details of each 
experiment, we point out a few general features here. 



• Firstly, the choice of initial parameters in each experiment avoids inter- 
section of the body momenta with fixed points. It is well known that the 
coadjoint orbits of the classical rigid body with distinct moments of inertia 
have saddle points at (0, ±7r, 0) (which are connected by four heteroclinic 
orbits) and centers at (±7r,0, 0) and (0,0, ±ir). The numerical solution 
does not become unstable, however, provided the time step is no larger 
than approximately 0.5. 

• The numerical round-off error in each representation depends upon the 
principle moments of inertia. This is shown after 10 time steps in figures 
7 and 8. For example, when I\ = I2 > I3, the error in the (i) spatial 
angular momentum is O(10 -8 ) and O(10~ n ) and (ii) energy is O(10 -7 ) and 
O(10 -10 ) for the respective body and spatial representations. When Ii > 
I2 > -^3) the same respective errors in the (i) spatial angular momentum 
are O(10- n ) and O(10- 14 ) and (ii) energy are O(10- 10 ) and O(10~ 13 ). 
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• The results comparing Dormand and Prince's explicit Rungc-Kutta method 
(0) (implemented in the odc45 routine) and the DMV algorithms should 
not be interpreted as a performance comparison. In each experiment, the 
odc45 solver was run at the smallest time step possible purely to provide a 
quantitative benchmark for the DMV algorithm. 

• We find a good agreement between the numerical results and the analytic 
solution for the rigid body and confirm conservation of spatial angular 
momentum and the Casimirs ||M||| and (det(I), \\I\\) for the body and 
spatial representations respectively. For the body representation of heavy 
top motion, the Casimir (M, T) is also correctly computed. Our numerical 
results qualitatively match those obtained by McLachlan and Zanna H^) 
and Cellcdoni and Safstrom (0) for the body representation of the rigid 
body and heavy top respectively. 

• We perform two coupled rigid body experiments in which two identical 
bodies are subject to the same initial conditions (i) but are initially at right 
angles to each other and (ii) are initially aligned with each other. In the first 
experiment, shown in figure 10 we observe non-periodic behaviour in the 
body angular momentum components caused by exchanges of momentum 
between the two bodies. The second component of the momentum changes 
the most, ranging from —1 to 0.8. In the second experiment, shown in 
figure 11, we recover a rigid body motion similar to that shown in figure 5, 
except that M 3 varies. 
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Figure 5. This figure compares numerical simulations of a rigid 
body over 1000 time steps for the case when the principal moments 
of inertia h = h > h (h = 2, I 2 = 2, 7 3 = 1). The top three 
graphs each show a component of the body angular momentum 
of rigid body motion. The bottom two graphs show the Casimirs 
\\M ||| and ||/||2, det(I) of rigid body flow in the respective body 
and spatial representations. The graph labelled (i) 'body DMV 
is the solution computed by the body DMV integrator, (ii) 'spa- 
tial DMV is the body frame translated solution computed by the 
spatial DMV integrator (which computes the angular momentum 
and moment of inertia in the spatial representation) (iii) 'ode45' 
is an explicit Runge-Kutta (4,5) integrated solution of the Eulcr- 
Arnold equations using the Matlab routine, ode45 and (iv) 'ana- 
lytic' is the analytical solution. The initial conditions for this simu- 
lation are the initial body angular momentum components given as 
Mi(0) =0.1, M 2 (0) = 0, M 3 (0) = 1. The top three graphs show 
that the DMV momentum matches the analytical solution which 
describes the rolling of a cone of constant angle in the body on a 
second cone of constant angle fixed in space {2lh . The 2 nd from 
bottom graph shows that the body DMV integrator precisely com- 
putes the Casimir ||M||| |2lh suggesting preservation of the rigid 
body Lie-Poisson structure and consequently that the DMV angu- 
lar momentum remains on the sphere. The bottom graph shows 
that the Casimirs ||/||2 and det(I) of the spatial DMV integrator 
are correctly computed suggesting that the Lie-Poisson structure 
on the dual of the semi-direct product lie algebra is also preserved. 
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Figure 6. This figure compares numerical simulations of a rigid 
body over 1000 time steps for the case when the principal mo- 
ments of inertia I\ > I 2 > h {h = 3.5, h = 2.5, I 3 = 2). The top 
three graphs each show a component of the body angular momen- 
tum of rigid body motion and the bottom two show the Casimirs 
\\M ||| and ||/||2, det(I) of rigid body flow in the respective body 
and spatial representations. The initial conditions for this simula- 
tion are the initial body angular momentum components given as 
Mi(0) = -0.5, M 2 (0) = 0, M 3 (0) = 1. The top three graphs show 
that the DMV momentum matches the analytical solution describ- 
ing the intersection of energy ellipsoids with coAdjoint orbits of 
SO(3) which are two-spheres |2lT) . Note that although our choice 
of simulation parameters avoids the flow intersecting either of the 
two saddle points at (0, ±| \M\ I2, 0) or centers at (±||M 1 12, 0, 0) and 
(0, 0, ±| \M\ I2), the solution does not become unstable provided the 
maximum time step is ~ 0.5. 
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Figure 7. This figure compares the energy and spatial angular 
momentum error in numerical simulations of rigid body motion 
over 10000 time steps for the case when the principal moments of 
inertia I\ = l<x > I3 {I\ = 2, I2 = 2, 23 = 1). The top graph shows 
the relative energy error growth in the solutions of the 'DMV and 
'ode45' integrators and the bottom graph shows the error in the 
approximated spatial angular momentum. The initial conditions 
for this simulation are the initial body angular momentum compo- 
nents given as M x (0) = 0.1, M 2 (0) = 0, M 3 (0) = 1. The graphs 
show that the error in the body DMV integrator computation of 
the energy and spatial angular momentum is higher than the error 
computed by the spatial DMV integrator. The bottom graph also 
shows that the spatial DMV integrator conserves spatial angular 
momentum to numerical round off. 
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Figure 8. This figure compares the energy and spatial angular 
momentum error in numerical simulations of rigid body motion 
over 10000 time steps for the case when the principal moments 
of inertia h > I 2 > h (h = 3.5, 7 2 = 2.5, / 3 = 2). The top 
graph shows the relative energy error growth in the solutions of 
the 'DMV and 'ode45' integrators and the bottom graph shows the 
error in the approximated spatial angular momentum. The initial 
conditions for this simulation are the initial body angular momen- 
tum components given as Mi(0) = -0.5, M 2 (0) = 0, M 3 (0) = 1. 
In contrast with the previous figure, the graphs show that the er- 
ror in the spatial DMV integrator computation of the energy and 
spatial angular momentum is higher than the error computed by 
the body DMV integrator. The bottom graph also shows that 
the body DMV integrator conserves spatial angular momentum to 
numerical round off. 
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Figure 9. This figure compares numerical simulations of the 
body representation of the heavy top over 1000 time steps for 
the case when the principal moments of inertia I\ = I2 > ^3 
(ii = 2, I2 = 2, I3 = 1). The top three graphs each show a compo- 
nent of the body angular momentum of heavy top motion and the 
bottom graph shows the error in the 'body DMV and 'ode45' com- 
puted Casimir (M, V) of the heavy top Lie-Poisson bracket. The 
initial conditions for this simulation are the initial (i) body angu- 
lar momentum components and (ii) position of the vertical axis 
in the body frame given respectively as Mi(0) = 0.1, M2(0) = 
0, M3(0) = 1 and T = [0,0,1]. Whenever the first and second 
components of the body angular momentum are non-zero, heavy 
top motion breaks the S 1 symmetry about the vertical axis. The 
bottom figure confirms that the Casimir (M, T) is always con- 
served. 
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Figure 10. This figure compares numerical simulations of the 
coupled rigid body, as seen in the frame of body 1, over 1000 
time steps for the case when the two identical rigid bodies are 
initially positioned at right angles to each other. The top three 
graphs each show a component of the body angular momentum of 
body 1 and 2. The initial conditions for this simulation are the 
initial (i) body angular momentum components (ii) orientation of 
the bodies relative to their E% axes and (iii) angle between the 
mechanical attachments at the ball and socket joint given respec- 
tively as M 2 (0) = M^O) = [0.5,0, 1], 0(0) = tp{0) and 0(0) = f . 
The principal moments of inertia of the two identical rigid bodies 
are I± = I2 > I3 (ii = 2, 1% = 2, I3 = 1). The graphs show 
that the components of body angular momentum are not periodic 
and extreme values are different from those of the single (uncou- 
pled) rigid body shown in figure 5, indicating transfer of angular 
momentum between the two bodies. 
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Figure 11. This figure compares numerical simulations of the 
coupled rigid body, as seen in the frame of body 1, over 1000 
time steps for the case when the two identical rigid bodies are ini- 
tially aligned with each other. The top three graphs each show 
a component of the body angular momentum of body 1 and 2. 
The initial conditions for this simulation are the initial (i) body 
angular momentum components (ii) orientation of the bodies rel- 
ative to their E% axes and (iii) angle between the mechanical 
attachments at the ball and socket joint given respectively as 
M 2 (0) = M^O) = [0.5,0,1], 0(0) = V(0) and 0(0) = 0. The 
principal moments of inertia of the two identical rigid bodies are 
h = h > h (h =2, I 2 = 2, 7 3 = 1). The graphs show that the 
components of body angular momentum of both bodies arc similar 
to those of the rigid body motion shown in figure 5, differing in 
the M3 component, which varies here. 
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Figure 12. This figure shows the error in computation of the 
Casimirs and conserved spatial angular momentum of the cou- 
pled rigid body motion, as viewed in the frame of body 1. for 
the case when both identical bodies are initially positioned at 
right angles to each other. The top graph shows the absolute 
error in 'DMV and 'ode45' computation of the C.R.B. Casimir 
\\M\\l = 1 1 Mi + AM 2 A T ||| The middle and bottom graphs 
show the comparative error in the energy and spatial angular mo- 
mentum of the coupled rigid body. The initial conditions for this 
simulation are the initial (i) body angular momentum components 
(ii) orientation of the bodies relative to their E3 axes and (iii) angle 
between the mechanical attachments at the ball and socket joint 
given respectively as M 2 (0) = M^O) = [0.5, 0, 1], 0(0) = ^(0) and 
0(0) = The principal moments of inertia of the two identical 
rigid bodies are I\~li> I3 (I\ = 2, I 2 = 2, I 3 = 1). 
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10. Conclusion 



The spatial and body representations ol rigid body motion correspond to the 
spatial and convective representations of continuum dynamics. The discrete Cleb- 
sch approach (0) provides a unified framework to derive variational integrators for 
both descriptions of the continuum. In this paper we demonstrate the utility of 
this framework by deriving explicit variational integrators for various rigid body 
motions in the spatial and body representations. We derive the momentum maps 
corresponding to the symmetry reductions of the discrete Lagrangians, the discrete 
EP equations and consequently prove, where appropriate, conservation of spatial 
angular momentum. 

This paper has pursued the relationship between variational integrators, derived 
in the discrete Clebsch framework, and existing studies of discrete integrable rigid 
body systems. In the body representation, we recover a discrete integrable analogue 
of the Euler- Arnold equations, first discovered by Moser and Veselov fl^l . In the 
spatial representation we obtain discrete EP equations with an adverted parameter. 
These correspond to the Lie-Poisson equations on the dual space of a semidirect 
product lie algebra discovered by Bobenko and Suris (0). This discovery provides 
a discrete extension to the work by Holm, Marsden and Ratiu l[l3h who developed 
the theory of EP entirely within a Lagrangian framework so that the EP equations 
with advection always correspond to Lie-Poisson Hamiltonian systems on the dual 
of a semi-direct product Lie-algebra. Consequently, the DMV equations have a 
family of Casimirs associated with the Lie-Poisson bracket for these systems, two 
of which we confirmed by numerical simulation of the spatial representation of the 
rigid body. 

We provide several numerical experiments to demonstrate the conservative prop- 
erties and accuracy of the explicit variational integrators that we derived in each 
case. We find a good agreement between the numerical results and the analytic 
solution for the rigid body and detect conservation of the spatial angular momen- 
tum and the total body angular momentum Casimir when the heavy top motion 
precesses purely about the vertical axis. Additionally, the numerical results in this 
paper qualitatively match those obtained by McLachlan and Zanna (|2fJ) for the 
rigid body and Celledoni and Safstrom (0) for the heavy top. 

We also observe non-periodic behaviour and a four-fold increase in the extremal 
values of the body angular momentum (relative to the uncoupled case) caused by 
exchanges of momentum between two identical coupled rigid bodies, which are both 
initialised with the same body angular momentum, but initially positioned at right 
angles to each other. 

We do, however, observe some less favourable and unexpected features in the nu- 
merical experiments. Firstly, in the rigid body experiments, the numerical solution 
becomes unstable and unconscrvative if the coAdjoint orbits on the sphere intersect 
either of the saddle points at (0, ±7r, 0) (which are connected by four heteroclinic 
orbits) or centers at (±7r, 0, 0) and (0, 0, ±7r) unless the time step is approximately 
no larger than 0.5. Precise bounds on the time step should be determined before 
applying DMV integrator to the study of bifurcating systems such as the heavy top 
by Lewis, Simo and Marsden |l8h and the dynamics and stability of coupled rigid 
bodies by Sreenath, Krishnaprasad and Marsden |29h and geometrically exact rod 
by Simo, Posbergh and Marsden (|28l) . 



Secondly, we find that the relative round-off error growth between the two rep- 
resentations is dependent upon the principle moments of inertia. For example, 
when I± = I2 > I3, the error in the (i) spatial angular momentum is O(10 -8 ) and 
O(10 -11 ) and (ii) energy is O(10~ 7 ) and O(10 -10 ) for the respective body and 
spatial representations. When I\ > I2 > I3, the same respective errors in the (i) 
spatial angular momentum are O(10 -11 ) and O(10~ 14 ) and (ii) energy are 0(1CP 10 ) 
and O(10 -13 ). For very long simulations, this aspect should not be overlooked and 
there may be practical benefit in choosing one representation over the other for the 
purpose of minimising round-off error. 

The discrete and continuous versions of the equations of motion and momen- 
tum maps are remarkably similar. Appendix A compares the body and spatial 
representations of these quantities in their continuous and discrete forms. 

In a forthcoming paper, we will attempt to address outstanding numerical issues 
and describe how the discrete Clebsch approach can be extended for ellipsoidal 
motions and motions of a geometrically exact rod in the convective representation. 
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Appendix A. Body and Spatial Representations in Continuous and 

Discrete Time 



Property 


Continuous 


Discrete 




Body attitude 


A(t) 


e so{N) 


A k e SO{N) 




Angular velocity 


n = 


a t a = -n T 


fifc+i = A^Afe + i 




Inertia Matrix 


Io 








Angular momentum 


M = 


= i n - n T i a 


M k = /ofife - nli 




Equations of motion 


M = 


= ad* n M 


M fe+ i = Ad* nT M k 




Right momentum map 


J R 


= PoA 


4 +1 = P k o A k 





Table 1. Comparison of the terms required to describe the body 
representation of the rigid body in continuous and discrete time 
as derived using the Clebsch approach. Blank items in the right- 
hand column indicate that they are identical to their discrete time 
descriptions. 





Property 


Continuous 


Discrete 






Body attitude 
Angular velocity 
Inertia Matrix 
Angular momentum 


A(t) € SO(N) 

LU — AA T — -LU T 

I = AJ A T 
m — Ilu — lu t I 


A fc e SO(N) 
ujk+i = Afc+iA^ 
Ik = A k IoA k 
m k = I k uj k - uj k I k 






Equations of motion 


m = ad^m — V/L o I = 0, 

i = [u,i] 


m k+ i = Ad* uT mk - Vi k LoI k , 

Ik+l = (^k+llk^k+l 






Left momentum map 


J L = PoA + Jol 


J L L +L =P k oA k + G k oI k 





Table 2. Comparison of the terms required to describe the spatial 
representation of the rigid body in continuous and discrete time. 
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Property 


Continuous 


Discrete 






Body attitude 
Angular velocity 
Inertia Matrix 
Angular momentum 
Orientation of the z-axis 


A(t) e SO(N) 
fi — A T A — -fi T 
Io 

M = iofi + fi T /o 

r = A T z 


A fc G SO(N) 
fife+i = Afc"Afc + i 

M k = /ofifc - fife Io 
T k = A£z 






Equations of motion 


M = ad^M + mgY o x, 

t = -fir 


Mfc+i = Ad* nT M k + mgTk. X, 
Tfc+i = fifc+iTfc 






Right momentum map 


j ii = PoA + roJ r 


Jk+i = Pk oA k +r fc o jfc 





Table 3. Summary of the terms required to describe the body 
representation of the heavy top in continuous and discrete time. 





Property 


Continuous 


Discrete 






Body attitude 
Angular velocity 
Inertia Matrix 
Angular mom. 
Position of c.o.m. 


A(t) G SO(N) 

CO = AA T = ~LO T 

I = A/ A T 
m = Ilj — w T I 
X = Axo 


Afc G SO(N) 
LOk+i = Afc+iAfc" 
h = A,, Jo A^ 
m k = huJk - ojklk 
Xk = AfcXo 






Eqns of motion 


rh — ad^m — ViL 1 + rngx o z, 

i = [u,r\, 

X = ux 


nik+i = Ad^ T m k - \/i k LoI k + mgxk o z, 

k ^ 

h+i = Wfc+i/fco; fc+1 , 
Xfc+i = Ufc+iXfc 






Left mom. map 


J L = AoP - J ol -3 x o x 


Jfc + i = Afc o Pfc + Gfc o I k +J^o X k 





Table 4. Summary of the terms required to describe the spatial 
representation of the heavy top in continuous and discrete time. 
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Property 


Continuous 


Discrete 






Body i attitude 
Body i angular velocity 
Position of c.o.m. of body i 
Orientat n matrix (rel. to b. 1) 

Mod. inertia matrix of body i 
Body 1 angular momentum 

Body 2 angular momentum 


Ai(t) G SO(N) 

n = Af A, = -nf 

di 

A \ T A 
A = A 1 A2 

J. j — Ii ' J-) i j 

Mi = AOi - flj h - eskew(Aft 2 -Di2) 
M 2 = I 2 fl 2 - fijf h + eskew(D 12 Qf A) 


A* S SO(N) 
0f +1 = Af A k+1 

Mf =i t QJi- olfh 

+2eskew(Di2 + Q k+lT ' D{ 2 ) 

m$ = m k 2 - nfi 2 

-2eskew(Di2 - n k+lT D 12 ) 






Equations of motion 


Mi = ad,Q.Mi, 
A = Afi 2 - fiiA 


M k+L =Ad* nkT M?, 
A fc+i = Q fe+i^ A fc fe+i 






Right momentum map 


J l R = Pi o Ai + J o A 


j£ = P k o A k + J k o A k 





Table 5 . Summary of the terms required to describe the body rep- 
resentation of the coupled rigid body in continuous and discrete 
time, where i G {1,2}. 





Property 


Continuous 


Discrete 






Body attitude 
Angular velocity 
Inertia matrix 
C.o.m. matrix 

Sp.a.mom.(Bl) 
Sp.a.mom.(B2) 


A(t) e SO(N) 
LU = AA T = -lu t 
h = AjfAj 
Dvt = A 1 D° 12 A 2 " 

mi = UJ1I1 — Iiut'i 
+ eske w ( D\ 2 oj 2 ) 

m 2 = lu 2 I 2 — I 2 u] 2 
— eskew(tJi_Di2) 


A fe G SO(N) 

Ldk+l = Afe+lA^ 

j|= = A k J°Af 
D k 2 = A\D\ 2 Af 

r fc+i 
1 1 

, A V 

m\ = 2skew{i k Lo k+1 - eto k+lT (uj k+1 - I d )D k 2 ) 

1 2 

, A V 

m\ = 2skew{I 2 Lo 2 +1 ~ eu,% +lT (u k+1 - h)D k 12 ) 






Eqns of motion 


m = Yli a ^i m i + V/.Lo/i, 
+X7 Dl2 LoDl 2 , 

Ii = [Ui, Ii], 

D12 = w\D\2 — D\ 2 uj 2 


m k+1 = T.rAdl kT m k+1 + V/*Lo J?, 

fk+l _ , fc+1 fk k+l T 
1 i - ^i 1 i w i 

U \2 — W l JJ 12 UJ 2 






Left mom. map 


J L =T, i Pi°A l + J l o It 
+JoD 21 


J L = £<■*?<> A- + G k o If + J fe o £> 2 fc i 





Table 6. Summary of the terms required to describe the spatial 
representation of the coupled rigid body in continuous and dis- 
crete time, where i G {1,2}. 
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Appendix B. The Spatial DMV Algorithm for the Rigid Body 

For completeness, we include the specification of the (explicit) DMV algorithm 
for the spatial representation of the rigid body. Following the approach of l|25r) . the 
algorithm uses the Schur decomposition of the Hamiltonian for the algebraic Ricatti 
equation to construct a symmetric matrix S k . Numerical experiments in sectional 
show that there is little difference between the conservative properties of the spatial 
and body versions of this algorithm. We observe that the numerical round-off error 
differs between the two different versions of the DMV algorithm depending upon the 
principal moments of inertia. Numerical experiments, not presented in this paper, 
find negligible difference in the stability and computational performance between 
the two versions. 

(1) For k = 1 :-> NT 

(2) m k = uj k _ 1 m k _ 1 ul_ 1 + [/*._!, w fe _i +u%_ x \ 

(3) Hk = ( (i' 2 , 4 ) 

(4) [R k , U k ] = Schur(?4) 

(5) S k = (R k hi(R k )n 

(6) u k = (S k + ^I, 1 

(7) h = LOkh-i^l 

(8) k = k + 1 
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